Symbolic computation of the Hartree-Fock energy from 
a chiral EFT three- nucleon interaction at N 2 LO 



q B. Gebremarianr a ' b , S. K. Bogner a,b , T. Duguet a ' b,c 

O 

"■National Superconducting Cyclotron Laboratory, 1 Cyclotron Laboratory, East-Lansing, 
q ML 48824, USA 

D b Department of Physics and Astronomy, Michigan State University, East Lansing, ML 

Q 48824, USA 

C CEA, Centre de Saclay, LRFU /Service de Physique Nucleaire, F-91191 Gif-sur-Yvette, 
_^ France 



Q_i Abstract 

We present the first of a two-part Mathematica notebook collection that 
q implements a symbolic approach for the application of the density matrix 

c/5 expansion (DME) to the Hartree-Fock (HF) energy from a chiral effective 

.^h field theory (EFT) three-nucleon interaction at N 2 LO. The final output from 

the notebooks is a Skyrme-like energy density functional that provides a 

quasi-local approximation to the nonlocal HF energy. In this paper, we 
£j discuss the derivation of the HF energy and its simplification in terms of 

the scalar /vector-isoscalar/isovector parts of the one-body density matrix. 

Furthermore, a set of steps is described and illustrated on how to extend the 

approach to other three-nucleon interactions. 
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Nature of problem: The calculation of the HF energy from the chiral EFT three- 
nucleon interaction at N 2 LO involves tremendous spin-isospin algebra. The prob- 
lem is compounded by the need to eventually obtain a quasi-local approxima- 
tion to the HF energy, which requires the HF energy to be expressed in terms of 
scalar/ vector- isoscalar/isovector parts of the one-body density matrix. The Math- 
ematica notebooks discussed in this paper solve the latter issue. 
Solution method: The HF energy from the chiral EFT three-nucleon interaction at 
N 2 LO is cast into a form suitable for an automatic simplification of the spin-isospin 
traces. Several Mathematica functions and symbolic manipulation techniques are 
used to obtain the result in terms of the scalar/vector-isoscalar/isovector parts of 
the one-body density matrix. 
Running time: Several hours. 

PACS: 02.30.GP; 02.60.Jh 

Keywords: three-nucleon interaction, Symbolic Hartree-Fock, Symbolic Density 
matrix expansion 

1. Introduction 

The nuclear Energy Density Functional (EDF), due to its computational 
advantages, is the only tractable approach to the calculation of ground- 
and excited-state properties of medium to heavy mass nuclei pp. Currently, 
nuclear EDF calculations rely on phenomenologically adjusted Skyrme and 
Gogny functionals. Even though these empirical functionals have been suc- 
cessful in providing a satisfactory description of bulk and certain spectro- 
scopic properties of known nuclei, they lack predictive power away from the 
region where they are fit. In light of this, various groups are pursuing dif- 
ferent strategies to improve the predictive power of such phenomenological 
functionals [2j El IH El El [7]. One of these approaches is based on the ex- 
plicit addition of microscopic long-range physics [H |9] through the use of 
the density matrix expansion method [TUl ITT] . 

Recent progress in evolving chiral effective field theory (EFT) interac- 
tions [121 [13] to lower resolution scale using the renormalization group 
methods [HI [151 HE] makes the construction of microscopically-based en- 
ergy density functionals a feasible endeavor. Towards such a goal, we de- 
rive the Hartree-Fock contribution from chiral EFT two- and three-nucleon 
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interaction at N 2 LO. Furthermore, the highly nonlocal nature of the result- 
ing energy density functional requires the application of the density matrix 
expansion [101 [XT] to derive a Skyrme-like quasi-local energy density func- 
tional [§1 [9] . Still, the EDF thus obtained only contains the Hartree-Fock 
physics such that further correlations must be added to produce any reason- 
able description of nuclei. In the short term, such an addition of correlations 
can be done empirically [TTj . Eventually though, it is our goal to compute 
higher-orders in perturbation theory and design a generalized density matrix 
expansion that is suited to simplify the non-localities in both space and time 
that arise beyond the Hartree-Fock level [18J. 

The tremendous algebra required to derive a Skyrme-like energy density 
functional from the exact HF energy through the DME makes a manual 
calculation impractical. This is especially true for the contribution from 
the three-nucleon interaction. In addition, such a derivation displays several 
features that make it amenable to symbolic automation [23]: (i) it involves 
many similar and repetitive algebraic steps (ii) most of it does not involve 
numerical computation, and (iii) the part of it that seems to require numerical 
computation, such as multidimensional integrals, can be performed using a 
combination of analytic expansion and symbolic integration [T9] . 

In this paper, we present the first part of the Mathematica solution in 
which we implemented the calculation of the HF energy from the chiral EFT 
three-nucleon interaction at N 2 LO and its subsequent simplification and re- 
expression in terms of scalar /vector-isoscalar/isovector parts of the one-body 
density matrix. The second part, which is the subject of a subsequent pub- 
lication, will deal with the application of the DME to the output of the first 
part, thereby yielding a Skyrme-like energy density functional. The paper 
is structured as follows: in section |2j the starting expression of the HF con- 
tribution from the chiral EFT three-nucleon interaction at N 2 LO is written 
in such a way that its Mathematica implementation becomes transparent. 
This is followed by section [3] where the Mathematica implementation is dis- 
cussed in detail. The subject of section [4] is devoted to the organization of 
the notebooks and the generation of the Mathematica code itself from the 
attached Python script, as well as to the modifications required to adapt the 
application of the notebook to more general three-nucleon interactions and 
to relax the isospin symmetry imposed in the present work. The last section 
contains the conclusions. 
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2. HF energy from the chiral EFT three- nucleon interaction at 
N 2 LO 

2.1. The \-EFT three-nucleon interaction at ISPLO 

From a general standpoint, three-nucleon interactions can be written as 

(k 1 k 2 k 3 |V A 3A r|k{k2k3) = — 6 kl+k2+k3 _ k /_ k ^_ k ^ V 3 jv(kik 2 k 3 |k(k2k 3 ) , (1) 

where Q is the volume used in the box-normalization of the momentum basis 
states |k), 5k 1 +k2+k 3 ,-kj'-k^-k ; ; is the Kronecker delta and ^^(kxk^lk^k^kg) 
is a matrix element in momentum space and an operator in spin-isospin space. 
In the above expression, the dependence of the interaction on spin and isospin 
degrees of freedom is not displayed. The chiral EFT three-nucleon interaction 
at N 2 LO has three main components (i) the E-term (ii) the D-term and (iii) 
the C-term that are diagrammatically represented in Fig. [T] As the present 
work relies only on the operator structure of the interaction which is provided 
below, we do not specify the actual values of the various constants that appear 
in the interaction. Refer to Ref. [20] for details. Furthermore, we neglect the 
regulator function as we are dealing with a HF calculation. This is justified 
because the HF calculation does not probe high-momentum single-particle 
states. Specifically, one can neglect the regulator so long as k F <C A, where 
kp is the local Fermi momentum and A is the cutoff parameter employed in 
the interaction. The absence of the regulator function makes the interaction 
local in coordinate space. Even though section [5] illustrates how to extend 
the approach to non-local interactions, the locality of the interaction brings 
significant reduction in the complexity of the second part of the problem, 
viz, the application of the DME to produce a quasi-local EDF. 

2.2. Remarks on the notation 

The following definitions hold for the various operators that appear in 
the analytical expressions for E-, D- and C-terms of the interaction that are 
stated below. er f = (er i)X , a i>y , a i}Z ) and n = (r i)X , T i>y , a ijZ ) refer to the spin 
and isospin Pauli matrices of the i th particle respectively. The momentum 
transfer for the i th particle is denoted with q« = k« — k'^ where kj and k'j are 
the out- and in-momentum coordinates of the particle. The scalar product 
of two vectors a and b is denoted with a • b. Furthermore, repeated indices 
imply summation over those indices (Einstein convention). 
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Figure 1: The two-pion exchange, one-pion exchange, and contact parts of chiral EFT 
three-nucleon interaction at N 2 LO and the relevant coupling constants of the interaction. 



2.2.1. The E-term 

The E-term, which is a three-nucleon contact interaction, is the simplest 
part of the interaction. The corresponding expression reads 



where 



V^(kik 2 k3|kik2 k 3) = E{ n • r 2 + r 2 • r 3 + r 3 ■ n ) , 

Ce 



E 



(2) 
(3) 



2.2.2. The D-term 

The D-term involves a one-pion exchange and a contact interaction. Its 
analytical form reads 



\/»(k;k 2 k,kfk.%) = ' /! CU 



ai ■ q 2 a 2 - Q2 

q% + ml 



Ti ■ t 2 + 



o"2 ■ q,3 o"3 • Q3 



cr 3 • qi <7i • qi 

x r 2 • r 3 H 2— - — r 3 ■ n 



(4) 



5. JTie C-term 

The C-term involves two-pion exchanges and its analytic form is 



Vc(kik 2 k3|ki^) 
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where 



KH = <W [ -4 ^ + 2 ^q, ■ q, ] + f 2 e^rl a k • (q, x q,) . (6) 

J TV J TV J TV 

2.3. Basic form of the HF energy from the \-EFT three-nucleon interaction 
at tfLO 

The three-nucleon interaction can in general be decomposed as a sum of 
three terms 

V 3N = Via + V 23 + V 13 , (7) 

where V^ is symmetric in i and j. Specifically, for the chiral EFT three- 
nucleon interaction at N 2 LO, Vy depends on qj and qj and, in general, on 
the spin-isospin coordinates of the three nucleons. After a few basic alge- 
braic manipulations, the HF energy from the three-nucleon interaction can 
be expressed in terms of only one of the three Vy operators, e.g. V23, as 

0^ F > = ^(ufc|V 2 3(l-2P 1 3-P 2 3 + 2P 12 P 23 )|uA;), (8) 

ijk 

where P\ m denotes the exchange operator (of particles I and m) whereas 
i,jandfc denote occupied HF single-particle states. Note that for ease of 
notation, we are using the single-particle basis that diagonalizes the one- 
body density matrix of the HF Slater determinant. P\ m is defined as 

Plm = Plm Plm P\m i (9) 

with P[ m the coordinate exchange operator and the spin-isospin exchange 
operators given by P° m = (1 + o x ■ a m )/2 and P{ m = (1 +ti ■ r m )/2. One can 
identify three groups of terms in Eq. ([8]): direct, single exchange and double 
exchange termd^J The direct term corresponds to the expectation value of 
V23, the single-exchange term to the expectation value of V^— 2P 13 — P 23 ) 
and the double-exchange term to that of 2 V23P12P23 



lr This should not be confused with one- and two-pion exchanges contribution to the 
three-nucleon interaction. 
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<v?/ lx ) 



l^(ijk\V2s(-2Pi3-P2s)\ijk) 



J2^ k \^3Pl2P2 3 \tjk). 

ijk 



(10) 

(11) 

(12) 



As the derivation of the a Skyrme-like quasi-local EDF from the exact HF 
energy requires the application of the DME, we need to express the HF energy 
in the |r) £§) \a) ® |r) single-particle basis. This is due to the fact that the 
DME, as formulated [TUl [TTj . is most naturally applicable to the one-body 
density matrix expressed in coordinate space. Hence, we perform inverse- 
Fourier transformation of the interaction to express the matrix elements in 
|r) (g> |c) ® |t) single-particle basis. This transformation leaves the spin-isospin 
dependencies untouched. Confining the discussion to the k dependence, we 
have 



(rir 2 r 3 |V23|r(r2 r 3) 



6(t! - r[) S(v 2 - r 2 ) 5(r 3 - r 3 ) 

x V 23 (r 2 - ri,r 3 - n) 



where 



(13) 



V 23 (r 2 -r 1 ,r 3 -r 1 ) = --^- J dq 2 dq 3 e **-fc»-i> ^-(rs-n) V 23 (q 2 , q 3 ) . (14) 

At this point, we do not actually perform the integrals over the momentum 
coordinates |^] in Eq (14). Rather, Eq. (14) is used used as it is, resulting 
in five-dimensional integrals in Eqs. ( 19 )-(21 ) that are performed after the 



application of the DME as discussed in Ref. [9] 



The next target is to rewrite Eqs. rtl0h-(12) in a form transparent for 



Mathematica implementation (in |r) ® \a) ® |r) single-particle basis). We 



2 Except for the E-tcrm of the interaction which is a trivial thrce-nucleon contact in- 



teraction, thereby yielding simple delta functions as shown in Eq. (27 1 
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illustrate the steps required to achieve that with Eq. (12), for which we have 



ijh 



r 3 3 

= II dT 'm n rfr «^ A; l r ' /l(j/ l r l r/2(T 2 r 2 r/ 34 r 3) 

ijk a' 1 ..a3 t[..tz m=l n=l 

x (r/icr^ y! 2 o' 2 t' 2 rf 3 cr' 3 r 3 \V 23 P[ 2 P^\ r i a i T i r 2 cr 2 r 2 r 3 a 3 r 3 ) 
x (ricrxn r 2 a 2 r 2 r 3 o 3 T 3 \P[ 2 P 23 \ijk) , (15) 

where we used completeness relations 

2j / EI dYi l 1 " 10 " 17 " 1 r 2°"2T 2 r 3 cr 3 T 3 ) (ricriT! r 2 cr 2 r 2 r 3 a 3 r 3 \ = 1 , (16) 



Cl ..<J3 71 ..T3 i=l 



and = Pi m Pi m - We split the particle exchange operator such that the 
coordinate part acts on the wave-functions while the spin-isospin piece is 
taken care of along with the interaction. Let Xj represent (rj(7jTj) such that 
the one-body density matrix reads as 

g(Xj, X fc ) = yir ; r)/; r r k a k r k ) = 2J ( r fc <M r j r i) > ( 17 ) 



where the sums is over occupied single-particle HF states. Making use of 
this, we define another quantity, which we call the auxiliary density matrix, 

as 

^(Xj-jXfc) = q (r^cr-r- , v^n) , (18) 

where z e {1, 2, 3}. Basically, the spin-isospin coordinates of this quantity are 
those of the i th particle. Applying the steps demonstrated in Eq. (15) and 
using Eqs. (18), (13)-(14), one can express the direct, single-exchange and 
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double-exchange parts of the three-nucleon interaction HF energy as 
1 



Tr!Tr 2 Tr 3 



(^5T 1X > = -Tr!Tr 2 Tr 3 



x V 23 (r 2 - ri,r 3 - 
J d ri dr 2 dr 3 q\X 3 , Xx) £? 2 (X 2 ) g 3 (X u X 3 



(19) 



xV 23 (r 2 - ri,r 3 - n) Pj 



i ITT 

13 



- - TnTr 2 Tr 3 



J dv l dv 2 dv z q\X x ) £? 2 (X 3 ,X 2 ) £? 3 (X 2 ,X 3 



xV 23 (r 2 - r x ,r 3 - r x ) R 



■UT 

23 



(20) 



(^ F ' 2X ) 



Tr x Tr 2 Tr 3 / di?idr 2 dr 



^ 1 (X 2 ,X 1 ) f? 2 (X 3 ,X 2 ) e 3 (X!,X 3 



xV 23 (r 2 -r 1 ,r 3 -r 1 )P 2 7P ] 
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(21) 



where ^(Xj) = £? 4 (Xj, X 3 -) and Tr^ refers to tracing over the spin and isospin 
coordinates of the i th particle. The key to understand the form of these 
equations is the splitting of the particle exchange operator, performed in 



Eqs. (15), that results in the spin-isospin coordinates of each particle to be 



grouped in a single auxiliary density matrix. These are the basic equations 
to be implemented directly in Mathematica. The next section shows that 
the implementation of these equations is transparent, which would not have 
been the case without the trick used to group the spin-isospin coordinates of 
each particle in a single auxiliary density matrix. 



3. Mathematica Implementation 



Starting from Eqs. ( 19 )-(21 ), the implementation consists of the following 



two components: (i) the automated tracing operation of spin-isospin matrices 
and (ii) the re-expression of the HF energy in terms of the scalar /vector- 
isoscalar/isovector parts of the one-body density matrix. 
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3.1. Automated Tracing 



The first task to be automated is the tracing operations in Eqs. (19) 



(21). This requires expressing the auxiliary density matrix in terms of its 
scalar /vector-isoscalar/isovector parts. We adopt the notation used in Ref. 
[2T] where for two vectors, v and w in isospin space, vow denotes their 



scalar product. Hence 

Q'ir'a'fl ro-iTi) = ^ (p l (r', r) 5 <CTi S^ Ti + pj(r', r) o t t[u 5^ 

+ s*(r', r) • a <ai 5 <n + s{{r\ r) • <j <ai o t <t )j .(22) 

with the scalar-isoscalar, scalar-isovector, vector-isoscalar and vector-isovector 
components are given, respectively, as 

pi(r',r) = ^ Q^r'aiTi, ro-iTi), (23) 

<JT 

p\(r',r) = ^ W/, rein) r TiT >, (24) 

cft't 

s o( r '> r ) = Yl ^(rViri.nr.Ti)^, (25) 
si(r',r) = ^(r'^WiTi)^^, (26) 



where <r aia i. = {ai\a\a'^) andr TiT ; = (Ti\r\T-), with the operators o = (a x , cr y , a z ) 
and r = (r x ,T y ,T z ). In contrast to the one-body density matrix, we, occa- 
sionally, refer to the scalar /vector-isoscalar/isovector components of the one- 
body density matrix as nonlocal densities. Note that we maintain the explicit 
reference to the index i as it will be useful in the Mathematica representation 
of the respective quantity. At this point, we identify all the basic quantities 
that need to be represented by their own Mathematica symbol. These are 
listed in table [Tj In this table, the index % plays the following two roles: 

(i) It identifies the coordinate dependence of the scalar /vector-isoscalar/isovector 
nonlocal densities. For instance, for the double-exchange part of the 
HF energy (as given in Eq.(|2~T|)), the three auxiliary density matrices 
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have the coordinate dependence 



^ 2 (Xj,X fc ) = £ 2 (X 3 ,X 2 ), 
^(X^Xfc) = ^ 3 (X 1 ,X 3 )., 

implying that the scalar/vector-isoscalar/isovector nonlocal densities 
derived from these will have their coordinate dependencies fixed ac- 
cordingly. For instance, the Mathematica symbol plO refers to the 
scalar-isoscalar nonlocal density extracted from the first auxiliary den- 



sity matrix. According to Eqs.(19)-(21), this entails that plO denotes 



the non/local density po( r i) for the case of direct contribution, po( r 3, r i) 
for the first part of the single-exchange contribution, po( r i) for the 
second-part of the single-exchange contribution and p ( r 2 ; r i) for the 
double-exchange contribution. 

(ii) One should be convinced that for the spin and isospin traces, the index 
i plays the role of particle-id. Hence, the spin and isospin matrices 
obtained from the auxiliary density matrices are to be grouped with 
the respective matrices obtained from the interaction and spin-isospin 
exchange operators. 

Once the symbols of the various quantities at play are identified, the next 
step is to define the composite symbols, viz, vectors and tensors. In the 
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following we list such definitions. 

In[l]: = iM = Identity Matrix [2] ; 

ln[2]: = ax = {{0, 1}, {1, 0}}; 

In[3]: = <ry = {{0, -i}, {i, 0}}; 

In[4]: = ^ = {{l, 0}, {0,-1}}; 

In[5]: = rx = {{0, 1}, {1, 0}}; 

In[6]: = ry = {{0, -i}, {i, 0}}; 

In[7]: = rz = {{l,0}, {0, -1}}; 

In [8]: — a — {ax, ay, az}; 

In [9]: = r = {tx, ry, tz}; 
In[10]: = q2 = {q2x, q2y, q2z}; 
In[ll]: = q3 = {q3x, q3y, q3z}; 
In[12]: = pll = {pllx, plly, p\\z)\ 
In[13]: = slO = {slOx, slOy, slOz}; 

In[14]: = sll = {{sllrrrr, sllxy, sllxz}, {sllyx, sllyy, sllyz}, 

{sllzx, sllzy, sllzz}}; 
In[15]: = p21 = {p21x, p21y, p21z}; 
In[16]: = s20 = {s20x, s20y, s20z}; 

In[17]: = s21 = {{s21xx, s2Lry, s21x^}, {s21yx, s21yy, s21yz}, 

{s21zx, s21zy, s21zz}}; 
In[18]: = p31 = {p31x, p31y, p31z}; 
In[19]: = s30 = {s30x, s30y, s30z}; 

In[20]: = s31 = {{s31xx, s31xy, s31xz}, {s31yx, s31yy, s31yz}, 
{s31zx, s31zy, s31zz}}; 

The actual tracing operation involves: (i) grouping together all spin- 
isospin operators that originate from the interaction, spin-isospin exchange 
operators or an auxiliary density matrix carrying the same particle-id. This 
should be done without mixing their actual ordering as only matrices with 
different particle-ids commute with each other, (ii) Performing the spin and 
isospin traces separately. In the following, we demonstrate the steps followed 
using the simplest term from the HF energy of chiral EFT three-nucleon 
interaction (at N 2 LO), viz, the direct part of the HF energy from the E-part 
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of the interaction. We denote this quantity by (V^ ,E ' dlv ) and, according to 
Eqs. @ and |L9]), it reads 



/T/ HF,£,dir\ 
\ V 3N I 



1 



ETnTr 2 Tr 3 / dr lC /r 2 dr 3 ^(Xx) £? 2 (X 2 ) g 3 (X 3 ) 



x 5(r 1 - r 2 ) 5(r 1 - r 3 ) r 2 • r 3 . 
The Mathematica code to perform the automated tracing reads 



(27) 



In[19]: 



V3NDirectHFE 



3 1 
^ 128 



al=l 



plO * Tr[iM] * Tr[iM] 

+Tr[iM] * Tr[(pll.r)] + Tr[(slO.<r)] * Tr[iM] 



3 3 



+ EE sll P]]P 2 i]^p]]]*% 

/31=1 /3 2 =1 

p20 * Tr[iM] * Tr[iM.r[[al]]] + Tr[iM] * Tr[(p21.r).r[[al]]] 
+Tr[(s20.o-)] * Tr[iM.r[[al]]] 

*Tr[r[[/34]].r[[al]]] 



3 3 



+ ^^s2ip3]][[/34]]Tr[a 

/33=1 /34=1 

p30 * Tr[iM] * Tr[iM.r[[al]]] + Tr[iM] * Tr[(p31.r).r[[al]]] 
+Tr[(s30.<r)] * Tr[iM.r[[al]]] 

3 3 

+ EE s31[[/35]][[/?6]]Tr[a[[/35]]] * Tr[rp6]].r[[al]]] 

/35=1 /36=1 



Assumptions — > AssumptionsList 



(28) 
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where the Mathematica symbol AssumptionsList contains 



In[19]: = AssumptionsList = {pl0==p20==p30, pllz==p21z==p31z, 
pll y ==p21y==p31y == 0, pllx==p21x==p31x == 0, 
S 10x==s20x==s30x, S 10y==s20y==s30y, 
sl0z==s20z==s30z, sllxx==s21xx==s31xx == 0, 
sllxy==s21xy==s31xy == 0, sllxz==s21xz==s31xz, 
sllyx==s21yx==s31yx == 0, sllyy==s21yy==s31yy == 0, 
sllyz==s21yz==s31yz, sllzx==s21zx==s31zx == 0, 
sllzy==s21zy==s31zy == 0, sllzz==s21zz==s31zz} . (29) 

While the tracing operation has been explained, the list enclosed in 
Assumptions — > warrants a few lines of explanation. The list contains two 
groups (i) a set of statements in which symbols are equated with each other. 
These are solely the symbols referring to various scalar /vector-isoscalar/isovector 
densities having the same coordinate dependence, which can be traced back 



to the occurrence of the two delta functions in Eq. (27). Note that this is due 
to the fact that the E-term is a three-nucleon contact interaction. For the D- 
and C-terms, a more complicated coordinate dependence occurs that results 
in there being no symbols equated in the respective AssumptionsList. (ii) 
The second group contains a list of symbols that are set to zero. It is based on 
the fact that the HF energy is isospin invariant, as is the starting interaction. 
Moreover, the present application forbids proton-neutron mixing, implying 
that single-particle wave-functions are eigenstates of the isospin-projection 
operator on the z axis. Hence, we enforce proton-neutron symmetry by per- 
forming rotation in isospin space Ref. [2T] using the conserved symmetry 
operator 

Upn = 2exp^-^7rf 3 ^ , (30) 

which aligns the quantities along the z axis in isospin space. Such a step 
yields a tremendous reduction in the size of the resulting expressions. Thus, 
in the final expression of the HF energy, only isoscalar and the z component 
of the isovector densities appear. 



Even though Eq. (28 ) corresponds to the simplest part of the problem, its 
Mathematica code requires a difficult mix of proper ordering of the operators 
and correct labeling of the dummy indices for the summation operations. 
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This becomes truly involved when dealing, for instance, with the double- 
exchange part of the HF energy from the C-term of the interaction. To 
address this problem, we have written a Python script that can be used to 
generate the required Mathematica code. Refer to section [4] for details. 

At this stage, the output from Mathematica is such that the various 
vectors and tensors occur in terms of their components which makes the 
expression too lengthy and impractical. To illustrate this particular point, 
we choose a small part of the single-exchange contribution from the D-term. 
Refer to section [4] for details on the organization of the notebooks. Even 
in this case, the output is much simpler than the typical outputs from the 
double-exchange of the D- and any of the C-terma^J Keeping that in mind, 
the output reads 

Out[19]: = ^plO ^ 2q3yq3z (-3s20zs30y - 3s20ys30z + s21zzs31yz 

+ s21yzs31zz) + 2q3x(q3y(-3s20ys30x - 3s20xs30y 
+ s21yzs31xz + s21xzs31yz) + q3z(— 3s20zs30x — 3s20xs30z 
+ s21zzs31xz + s21xzs31zz)) + q3x 2 (-3s20xs30x + 3s20ys30y 
+ 3s20zs30z + s21xzs31xz — s21yzs31yz — s21zzs31zz — 3p20p30 
+ p21zp31z) + q3y 2 (3s20xs30x - 3s20ys30y + 3s20zs30z 

— s21xzs31xz + s21yzs31yz — s21zzs31zz — 3p20p30 + p21zp31z) 
+ q3z 2 (3s20xs30x + 3s20ys30y - 3s20zs30z - s21xzs31xz 

— s21yzs31yz + s21zzs31zz — 3p20p30 + p21zp31z" 
1 



-pll I 2q3yq3z (— 3s21xzs30y — 3s21yzs30z + s20zzs31yz 
8 V 

+ s20ys31zz) + 2q3x(q3y(— 3s21yzs30x — 3s21xzs30y 

+ s20ys31xz + s20xs31yz) + q3z(— 3s21xzs30x — 3s21xzs30z 

+ s20zs31xz + s20zs31zz)) +...]. (31) 



3r 



i The output of the automated tracing from the double-exchange of the D-term is tens 
of pages long while that of the single- and double-exchange of the the C-term are much 
longer than that. 
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Hence, the manual simplification of such an expressions is not feasible. The 
next section describes the technique we used to automate the task of writing 
these expressions in terms of the scalar /vector-isoscalar/isovectoi[^] parts of 
the density matrix: {p , p 1; s , Sx}. 

3.2. Scalar /vector-isoscalar /isovector re-expression 

The main objective of this part of the code is to rewrite the output of 
the automated parsing step in terms of the scalar/vector- isoscalar /isovector 
parts of the one-body density matrix with the proper coordinate dependence. 
The approach consists of two basic steps. Initially we form the eight groups 
that arise from choosing the scalar or vector components of for each of the 
three density matrices that are multiplied, which we represent schematically 
as 

groups = {{0|1}®{0|1}®{0|1}}, (32) 

where {0|1} represents taking one of the scalars {p , Pi} or one of the vectors 
{s , Si} of that density matrix. Once a particular group is formed, the next 
step involves forming the list of all scalars and vectors including momentum 
vectors coming from the interaction. 

In the following, a complete illustration is given for the double-exchange 
of the D-term. The list of scalars and vectors that we have for instance, for 
the {011} group of the double-exchange of the D-term, are 

groupOll = {q 3 , q 3 , pj, p\, sjj, s?, sjj, sf , } , (33) 

where the superscripts {1,2,3} specify the id of the original density ma- 
trix. Note that, for each density matrix, we have two possible choices 
(isoscalar /isovector) even after the scalar/vector choice has been made. In 
this set, the only momentum vector that appears is q% and it appears twice 
as can be seen from the combination of Eqs. Q and (21). 



Since energy is invariant with respect to rotation in coordinate, spin and 
isospin spaces, we next form the list of all possible scalars from the selected 
group. In each scalar, we can pick only one of the two possible choices orig- 
inating from the same density matrix. For instance, referring to groupOll, 



4 From here onwards, we use the term isovector freely even though we have set the 
first and second components to zero (only the third is nonzero). Refer to the discussion 
immediately following Eq. ( 29 I . 



16 



we cannot have pj and p\ in the same scalar. Hence, the scalars that can be 
formed from groupOll are 

scalarsOll = {pi (q 3 .q 3 ) (sQ.sg), pi (q 3 .q 3 ) (4-4), p} (q 3 .q 3 ) (sg.s?), 
pi (q 3 -q 3 ) (s?.sjj), pj (q 3 .sg) (q 3 .sg), p* (q 3 .s?) (q 3 .s?), 
pi(q 3 .sS)(q 3 .a?) ) Pi (03-s?) (qs-sg) }• (34) 

The subsequent step involves using Mathematica's SolveAlways function 
to obtain the coefficients for each term of scalarsOll. To complete the dis- 
cussion, let us explicit the part of the Mathematica code that implements 
the key ideas of this section 

in[19]:sllz = {sllxz, sllyz, sllzz}; 
in[19]:a21z = {s21xz, s21yz, s21zz}; 
in[19]:s3l£ = {s31xz, s31yz, s31zz}; 

in[19] :V3NDoubleExchangeHFDgroup011 = Simplify[V3NDoubleExchangeHFD, 
Assumptions — > {slOx == 0, slOy == 0, slOz == 0, 
sllxz == 0, sllyz == 0, sllzz == 0, 

p20 == 0, p21z == 0, p30 == 0, p31z == 0}] , (35) 

where V3NDoubleExchangeHFD is the symbol for the double-exchange con- 
tribution from the D-term obtained after the automated tracing. In the ac- 
tual implementation, the sheer size of the V3NDoubleExchangeHFD forces 
one to break the expression further into subcomponents. Refer to section 
[1] for details. In the above code, the definitions for {sllxz, sllyz, sllzz} 
allows one to collect the nonzero parts of the respective tensors, viz, the z 



components in isospin space as discussed in section 3.1 in to a vector. The 
code in Assumptions — > block simply sets the parts of the density matrices 
that are zero in the chosen group, viz, groupOll. For details, refer to the 
submitted Mathematica code. Once we obtain the particular group of inter- 
est, in this case groupOll, the coefficients of the various scalar terms listed 



in Eq. (34) are calculated using the code 
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in[19] :CoefficientList011 = SolveAlways \V3NDoubleExchangeHFDgroup01 1 

- (alg3.g3pl0s20.s30 + a2g3.g3pl0s21z.s31z + a3g3.g3pllzs20.s31z 
+ a4g3.g3pllzs21z.s30 + a5pl0g3.s20g3.s30 + a6pl0g3.s21zg3.s31z 
+ a7pllzg3.s20g3.s31z + a8pllzg3.s21zg3.s30) ==0, 
{plO, p20, p30, pllz, p21z, p31z, g3x, q3y, g3z, slOx, slOy, 
slOz, s20x, s20y, s20z, s30x, s30y, s30z, sllxz, sllyz, 
sllzz, s21xz, s21yz, s21zz, s31xz, s31yz, s31zz}] (36) 

The output symbol CoefficientListOll will contain the the values of the co- 
efficients {ai, a2, ....ag}- In this case 



in[19]:CoefficientList011 ={al -3/16, a2 1/16, a3 -> -3/16, 

a4 -> 1/16, a5 -> 3/8, a6 -> -1/8, 
a7^ 3/8, a8 -> -1/8}. (37) 

Finally a part of the code does the automatic replacement of the coordinate 
dependencies. As for the example used in the present section, one is dealing 
with a part of the double-exchange from the D-term, which implies that the 



coordinate dependencies of the three density matrices, according to Eq. (21 ), 
are 

^(X^Xfc) = g 1 (X 2 ,X 1 ), 
a 2 (X,,X fc ) = a 2 (X 3 ,X 2 ), 
a 3 (X J; X fc ) = ^(Xj.Xs). 

Thus far the multidimensional integrals and the interaction constants that 



we have in Eqs.(19)-(21 ) did not take part in any of the symbolic manip- 
ulations. Once the proper coordinate dependencies are restored, we insert 
the interaction constants and the multidimensional integrals resulting in the 
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final expression for V3NDoubleExchangeHFDgroup011, denoted as (Ion) 

P 7 

-3po(r2,riX(r3,r 2 )so(r 2 ,r 3 ) + p (r 2 , r a )af(r 3 , r 2 )4(r 2 , r 3 



/T/D2x\ _ -PA Cp 1 / , , /* 1 , p iq3.(r 3 -r 2 ) ^3 ^3 



- 3pi(r 2 ,ri)s£(r 3 ,r 2 )4(r 2 ,r 3 ) + pi(r 2 , ri)s^(r 3 , r 2 )so(r 2 , r 3 ) 

+ 3p (r 2 ,r 1 )s 3 (r 3 ,r 2 )s2(r 2 ,r 3 ) - p (r 2 , r^sf (r 3 , r 2 )s?(r 2 , r 3 ) 

+ 3pi(r 2 ,r 1 )s 3 (r3,r 2 )s]'(r 2 ,r 3 ) - px(r 2 , n)sf (r 3 , r 2 )s^(r 2 , r 3 ) 

(38) 

To find the total value of the contribution from the double-exchange of the 
D-term, similar analysis is carried out for the remaining seven groups and 



their results added to Eq. (38). Naively counting, one has to calculate about 
seventy-two group^Jof varying sizes to get the total contribution of the three- 
nucleon interaction. This shows the tremendous size of the algebra required 
to calculate the total HF energy and hence the need for automation. 

4. More on the Mathematica notebooks 

4-1. Automatic generation of the Mathematica codes and organization of the 
notebooks 



As mentioned in section 3.1 , writing the Mathematica codes is an error- 
prone procedure. This is mainly due to the fact that one has to include 
a large number of summation indices in a single expression and maintain 
proper operator orders. Consequently we wrote a python script that can 
be used to generate the Mathematic codes automatically given the proper 
parameters. Refer to the submitted files for details. 

The Mathematica implementation of the calculation of the HF energy 
and re-expression in terms of scalar /vector-isoscalar/isovector parts of the 
density matrix is carried out in eight notebooks. The contribution from the 



For the three parts of the interaction, we have direct, single-exchange and double- 
exchange with each requiring eight groups. 
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E- is implemented in a single notebook and the same for the D-term. The C- 
term is implemented in the remaining six notebooks. Refer to the submitted 
code for details. 



5. Possible Extensions 



In this section, we discuss the steps that can be taken to extend the 
approach so that one can treat non-local interactions and allow for proton- 
neutron mixing. 

5.1. Treating non-local interactions 

The approach can be extended to treat nonlocal interactions. The starting 



Eqs.( 19 )-(21 ) need to be modified for nonlocal interactions as 

3 3 



<v?/' dir > 



TriTr 2 Tr 3 



<v 3 H /' lx > 



TriTr 2 Tr 3 



1 



TriTr 2 Tr 3 



J] dr t H dv\ g 1 (X 1 , X[) g 2 (X 2l X' 2 ) £? 3 (X 3 , X' 3 ) 
i=i j=i 

x(rir 2 r 3 |V 23 |rir' 2 i/3) , (39) 

3 3 

d J] dv % \{ dv\ g\X 3 , Xi) £? 2 (X 2 , X' 2 ) g 3 (X 1 , X' 3 ) 

i=l j=l 



x(r 1 r 2 r 3 |V 23 |rir' 2 i/ 3 )P 1 7 



J] dv t J] dx\ g l {X u Xi) £? 2 (X 3 , X' 2 ) £? 3 (X 2 , X' 3 

i=l 3=1 



(v 3 T 2x ) 



TriTr 2 Tr 3 



3 3 

n ^ n dr i ^( x2 ' x i) ^ 2 ( x 3' ^ 3 ( x i' x s. 

i=l j=l 



(40) 



(41) 



Starting with this, the same automated spin-isospin tracing, scalar /vector- 



isoscalar/isovector re-expression can be carried out as described in sees. 3.1 



and |3.2| This easy extension remains valid throughout the first part of the 
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problem, specifically, calculation and re-expression of the HF energy in terms 
of scalar /vector-isoscalar/isovector parts of one-body density matrix. This 
is in contrast to the second part of the problem, application of the DME to 
produce a Skyrme-like energy density functional, where the non-locality of 
the interaction makes the problem significantly harder. 

The traditional Skyrme phenomenological interactions [22] employed in 
nuclear structure calculations are much easier to treat as they are zero-range 
non-local interactions. Hence, we demonstrate how the present symbolic 
approach can be applied directly to calculate the energy density functional 
from the three-nucleon part of Skyrme interactions. From the outset, it 
should be clear that we don't need to further apply the DME to obtain 
a local energy density functional as the interaction is of zero-range. The 
situation is partly similar to the E-term of the realistic chiral EFT three- 
nucleon interaction. However, the fact that we have gradient operators in 
the traditional Skyrme phenomenological interactions make the derivation 
harder than the corresponding derivation for the E-term. 

We start with a schematic three-nucleon Skyrme-like interaction 

(rir 2 r 3 |V 23 | r^r^) = (k£g + k' 2 | ) <J(n - ri)<J(r 2 - r' 2 ) 6(r 3 - r' 3 ) 

x 5(n - r 2 ) 6{r 2 - r 3 ) , (42) 

where according to the usual notation of Skyrme interactions = i/2(V< — 
Vj) and k^- = —i/2(V i — VQ. Our convention is that, when using such 



an interaction in Eqs. (39)- (41), gradients must act on the density matrices 
prior to making use of the delta functions. 

The next step in the symbolic simplification is to form the set of rules 
based on the relationship between the nonlocal scalar /vector-isoscalar/isovector 
densities with the various local densities that arise in Skyrme energy density 
functionals. A list and properties of these local densities and the relationship 
between these local densities and the nonlocal one they originate from can 
be found in Ref. [21]. For instance 

= V 2 po(R) - 2r (R), (43) 



3' 



=R 



where po(R) is the local scalar-isoscalar density and r (R) is the local scalar- 
isoscalar kinetic density. These relationships can be expressed as a Math- 
ematica replacement rules involving the symbols representing these quanti- 
ties. What follows is a chain of replacement rule operations to obtain the 
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final form of the local energy density functional. To illustrate the complete 
set of steps, we have coded and submitted a Mathematica notebook that 
derives the local energy density functional resulting from the the direct part 



(Eq. (39)) of the schematic three- nucleon Skyrme-like interaction given in 
Eq.(42). The most important advantage of this automation is that once the 
rules for mapping the nonlocal densities onto the local ones are established, 
calculating the energy density functional for different spin-isospin and gradi- 
ent structures requires a negligible fraction of the effort required to perform 
the corresponding derivation manually. 

5.2. Proton-neutron mixing 

Allowing for proton-neutron mixing makes possible to produce a more 
general skyrme-like energy density functional. The merit of the present ap- 
proach is that this can be done directly by removing the second group of 



assumptions discussed just after Eq. (29). The main negative side effect 
that one faces in doing so relates to the size of the resulting expressions 
that increase significantly, resulting in a possible failure of Mathematica's 
SolveAlways method. Usually, this problem can be alleviated by breaking 
the unhandled expressions into smaller pieces that can always be solved by 
SolveAlways. Moreover, allowing for proton-neutron mixing forces substan- 
tial changes when re-expressing the formula in terms of the scalar /vector- 
isoscalar/isovector parts of the one-body density matrix. The reason is that 
one has to use the full scalar-isovector, pij, and vector- isovector, s^-, symbols 
instead of their projections along the z axis in isospin space. As a result, one 
has to take this into account when constructing all possible scalars as de- 



scribed in section 3.2 Another approach can be to perform the calculations 
with a conserved proton-neutron symmetry, as is done in the present work, 
and map the resulting expression to the case of proton-neutron mixing. 



6. Comments and conclusions 

The first of a two-part symbolic approach to the derivation of a Skyrme- 
like quasi-local energy density functional from the HF energy of the chiral 
EFT three-nucleon interaction at N 2 LO was discussed and implemented in 
a collection of Mathematica notebooks. We have demonstrated that the 
tremendous spin-isospin algebra required to express the HF energy in terms of 
the coordinate-space scalar /vector-isoscalar /isovector parts of the one-body 
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density matrix can be effectively tackled by using Mathematica. In the sec- 
ond part of the work, which will be the subject of a subsequent publication, 
we show that the further application of the density matrix expansion method 
needed to generate the Skyrme-like energy density functional can also be 
completely automated. Possible future directions of research could be: (i) as 
a whole, the current Mathematica implementation conserves proton-neutron 
symmetry. Hence, one can allow proton-neutron mixing by relaxing this con- 
straint, thereby arriving at a more general Skyrme-like energy density func- 
tional. The steps that are required to achieve this have been outlined in the 
present contribution, (ii) We have demonstrated that the approach can be 
generalized to treat nonlocal interactions such as the quasi-local Skyrme in- 
teractions. Calculating the energy density functional starting from a Skyrme- 
interaction (including two- and three-nucleon parts) with a rich spin-isospin 
and gradient structures is a tedious task. In light of the recent interest along 
these lines (see Ref. [25J and references therein), extending the scheme to treat 
this problem could be useful, (iii) One can envision expanding the present 
work in such a way that first-order pairing correlations are treated along with 
the HF part, viz, performing HFB (Hartree-Fock-Bogoliubov) calculations. 
Combining this extension with proton-neutron mixing, one can have a start- 
ing Skyrme-like functional that can be used to handle proton-neutron pairing 
correlations as discussed in Ref. [21] • (iv) Implementing a similar scheme to 
treat four-nucleon interactions can also be one area of extension. 

7. Acknowledgements 

We thank R. J. Furnstahl for useful discussions. This work was supported 
in part by the National Science Foundation under Grant No. PHY-0758125 
and the UNEDF SciDAC Collaboration (DOE Grant DE-FC02-07ER41457). 

References 

[1] M. Bender, P. H. Heenen and P. G. Reinhard. Rev. Mod. Phys. 75 
(2003), pp. 121. 

[2] T. Lesinski, K. Bennaceur, T. Duguet and J. Meyer. Phys. Rev. C74 
(2006), pp. 044315. 

[3] T. Lesinski, M. Bender, K. Bennaceur, T. Duguet and J. Meyer. Phys. 
Rev. C76 (2007), pp. 014312. 



23 



[4] J. Margueron, H. Sagawa and K. Hagino. Phys. Rev. C77 (2008), pp. 
054309. 



[5] 
[6] 

[7 



[9 
[10 

[11 

[12 
[13 

[14 

[15 
[16 
[17 
[18 
[19 



T. Niksic, D. Vretenar and P. Ring. Phys. Rev. C78 (2008), pp. 034318. 

B. G. Carlsson, J. Dobaczewski and M. Kortelainen. Phys. Rev. C78 
(2009), pp. 044326. 

S. Goriely, S. Hilaire, M. Girod and S. Peru. Phys. Rev. Lett. 102 (2009), 
pp. 242501. 

S. K. Bogner, B. Gebremariam and T. Duguet. In preparation. 
B. Gebremariam, S. K. Bogner and T. Duguet. In preparation. 
J. Negele and D. Vautherin. Phys. Rev. C5 (1972), pp. 1472. 



B. Gebremariam, T. Duguet and S. K. Bogner. |arXiv:0910.4 979vl [nucl- 
th]. 

D. R. Entem and R. Machleidt. Phys. Rev. C68 (2003), pp. 041001(R). 

E. Epelbaum, W. Glockle and U. G. MeiBner. Nucl. Phys. A747 (2005), 
362. 

S. K. Bogner, T. T. S. Kuo and A. Schwenk. Phys. Rept. 386, 1 (2003); 



S. K. Bogner, A. Schwenk, T. T. S. Kuo and G. E. Brown, [ArXiv:nucT 
|th/0111042 vl. 

S. K. Bogner, R. J. Furnstahl, S. Ramanan and A. Schwenk. Nucl. Phys. 
A784 (2007), 79. 

S. K. Bogner, R. J. Furnstahl, S. Ramanan and A. Schwenk. Nucl. Phys. 
A773 (2006), 203. 

M. Stoitsov, M. Kortelainen, W. Nazarewicz, T. Lesinski, N. Schunck, 
S. Bogner, T. Duguet, D. Furnstahl and B. Gebremariam. Unpublished. 

V. Rotival, B. Gebremariam, T. Duguet, S. K. Bogner and R. J. Furn- 
stahl. Unpublished. 



B. Gebremariam, S. K. Bogner and T. Duguet. arXiv:0910.4993vl 
[physics, comp-ph] . 

24 



[20] E. Epelbaum. [arXiviOSl 1 . 1 133bvl [nucl-th] . 

[21] Perliriska, E. and Rohozinski, S. G. and Dobaczewski, J. and Nazarewicz, 
W. Phys. Rev. C69, 014316 (2004). 

[22] T. H. Skyrme. Phil. Mag. 1, 1043 (1956). 

[23] S. Wolfram. The Mathematica Book. Wolfram Media/ Cambridge Uni- 
versity Press, Cambridge, 1996. 

[24] Y. M. Engel, D. M. Brink, K. Goeke, S. J. Krieger and D. Vautherin. 
Nucl. Phys. A249 (1975), pp. 215. 

[25] J. Margueron and H. Sagawa. J. Phys. G: Nucl. Part. Phys. 36 (2009) 
125102. 



25 



Table 1: List of basic quantities and their corresponding symbols. 
Quantity Symbol Remark 



a 



Sjj sij 



The vector containing the three Pauli matrices: 
a = {cr x ,a y ,a z }. 

The vector containing the three Pauli matrices: 

The vector containing the components of the vector 
density: s = {s x , s y , s z }. The index ie{l,2, 3} denotes 



the first, second or third density in Eqs.(19)-(21) while 
j e {0, 1} refers to the isoscalar and isovector components 
respectively. For example, the Mathematica symbol slO 
defines the vector slO = {slOx, slOy, slOz}. 

The scalar part of the density matrix. The index 
ie{l, 2, 3} denotes the first, second or third 



Pij pij density in Eqs.(19)-(21) while je{0, 1} refers to 



the isoscalar and isovector parts, respectively. 
For example, the Mathematica symbol pll 
defines the vector pll = {pllx, plly, pllz}. 

The momentum vector q vector in the interaction 
qi where i e {2, 3} labels the particle. For example, 

the symbol q2 defines the vector q2 = {q2x, q2y, q2z}. 
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